"""
# Allen-Cahn Equation - Phase Field Modeling

This module implements the Allen-Cahn equation, a fundamental model in phase field 
theory for describing interface dynamics, phase separation, and ordering phenomena.
It's widely used in materials science, image processing, and mean curvature flow.

## Mathematical Background

The Allen-Cahn equation is:
```
∂φ/∂t = -M[δF/δφ] = M[ε²∇²φ - f'(φ)]
```

Where:
- φ: Order parameter (phase field)
- M: Mobility coefficient
- ε: Interface width parameter
- f(φ): Free energy density (typically double-well potential)

## Standard Double-Well Potential
```
f(φ) = (φ² - 1)²/4
f'(φ) = φ(φ² - 1) = φ³ - φ
```

## Key Features

- Interface motion by mean curvature
- Energy dissipative dynamics
- Conservation of total energy (not mass)
- Sharp interface limit as ε → 0

## Applications

- Phase separation dynamics
- Grain boundary motion
- Interface tracking
- Image denoising and segmentation
- Minimal surface problems

## References

- Allen, S. M. & Cahn, J. W. (1979). Acta Metall. 27, 1085.
- Evans, L. C. et al. (1992). Calc. Var. 1, 253-276.
- Fife, P. C. (1988). Dynamics of Internal Layers and Diffusive Interfaces.
"""

using LinearAlgebra, FFTW, Plots, Statistics, Printf
using HDF5  # For data export

export AllenCahnSystem, run_allen_cahn, analyze_interfaces, track_energy_dissipation

"""
    AllenCahnSystem

Configuration for Allen-Cahn equation simulations.

# Fields
- `mobility::Float64`: Mobility coefficient M
- `epsilon::Float64`: Interface width parameter ε
- `potential::Symbol`: Type of potential (:double_well, :obstacle, :custom)
- `potential_params::NamedTuple`: Parameters for potential function
- `domain_size::Tuple`: Spatial domain dimensions
- `grid_size::Tuple`: Grid resolution
- `boundary_condition::Symbol`: Boundary conditions (:periodic, :neumann)
"""
struct AllenCahnSystem
    mobility::Float64
    epsilon::Float64
    potential::Symbol
    potential_params::NamedTuple
    domain_size::Tuple{Float64, Float64}
    grid_size::Tuple{Int, Int}
    boundary_condition::Symbol
end

"""
    AllenCahnSystem(; mobility=1.0, epsilon=0.1, potential=:double_well, 
                    domain_size=(2.0, 2.0), grid_size=(128, 128), 
                    boundary_condition=:periodic, kwargs...)

Create Allen-Cahn system configuration.
"""
function AllenCahnSystem(; mobility=1.0, epsilon=0.1, potential=:double_well,
                        domain_size=(2.0, 2.0), grid_size=(128, 128),
                        boundary_condition=:periodic, kwargs...)
    
    # Default parameters for different potentials
    if potential == :double_well
        default_params = (a=1.0, b=1.0)
    elseif potential == :obstacle
        default_params = (barrier_height=1.0,)
    elseif potential == :custom
        default_params = NamedTuple(kwargs)
    else
        throw(ArgumentError("Unknown potential: $potential"))
    end
    
    potential_params = merge(default_params, NamedTuple(kwargs))
    
    return AllenCahnSystem(mobility, epsilon, potential, potential_params,
                          domain_size, grid_size, boundary_condition)
end

"""
    potential_derivative(φ, system::AllenCahnSystem)

Compute derivative of the free energy potential f'(φ).
"""
function potential_derivative(φ, system::AllenCahnSystem)
    if system.potential == :double_well
        # f(φ) = (φ² - a²)²/(4b²), f'(φ) = φ(φ² - a²)/b²
        p = system.potential_params
        return φ .* (φ.^2 .- p.a^2) / p.b^2
        
    elseif system.potential == :obstacle
        # Obstacle potential: f'(φ) = 0 for |φ| ≤ 1, ∞ otherwise (projected)
        return clamp.(φ, -1.0, 1.0) .- φ
        
    else
        throw(ArgumentError("Potential derivative not implemented for $(system.potential)"))
    end
end

"""
    potential_energy(φ, system::AllenCahnSystem)

Compute the free energy potential f(φ).
"""
function potential_energy(φ, system::AllenCahnSystem)
    if system.potential == :double_well
        p = system.potential_params
        return (φ.^2 .- p.a^2).^2 ./ (4 * p.b^2)
    elseif system.potential == :obstacle
        # Indicator function for |φ| ≤ 1
        return ifelse.(abs.(φ) .<= 1.0, 0.0, Inf)
    end
end

"""
    initialize_allen_cahn(system::AllenCahnSystem, initial_type::Symbol=:random)

Initialize the phase field with appropriate initial conditions.

# Arguments
- `initial_type::Symbol`: Type of initial condition
  - `:random`: Random field between -1 and 1
  - `:two_phase`: Two distinct phases
  - `:circular_interface`: Circular interface
  - `:sinusoidal`: Sinusoidal perturbation
  - `:droplets`: Multiple circular droplets
"""
function initialize_allen_cahn(system::AllenCahnSystem, initial_type::Symbol=:random)
    nx, ny = system.grid_size
    Lx, Ly = system.domain_size
    
    # Create coordinate grids
    x = range(0, Lx, length=nx)
    y = range(0, Ly, length=ny)
    X = [xi for xi in x, yi in y]
    Y = [yi for xi in x, yi in y]
    
    if initial_type == :random
        # Random initial condition
        φ = 0.1 * (2 * rand(nx, ny) .- 1)
        
    elseif initial_type == :two_phase
        # Sharp interface between two phases
        φ = ifelse.(X .< Lx/2, -0.9, 0.9)
        # Add small perturbation
        φ .+= 0.1 * sin.(2π * Y / Ly) .* exp.(-((X .- Lx/2) / (Lx/10)).^2)
        
    elseif initial_type == :circular_interface
        # Circular interface
        center_x, center_y = Lx/2, Ly/2
        radius = min(Lx, Ly) / 4
        distance = sqrt.((X .- center_x).^2 + (Y .- center_y).^2)
        φ = tanh.((radius .- distance) / (2 * system.epsilon))
        
    elseif initial_type == :sinusoidal
        # Sinusoidal perturbation
        φ = 0.1 * sin.(2π * X / Lx) .* cos.(2π * Y / Ly)
        
    elseif initial_type == :droplets
        # Multiple droplets
        φ = -0.9 * ones(nx, ny)
        n_droplets = 5
        for _ in 1:n_droplets
            center_x = Lx * rand()
            center_y = Ly * rand()
            radius = (min(Lx, Ly) / 10) * (0.5 + 0.5 * rand())
            distance = sqrt.((X .- center_x).^2 + (Y .- center_y).^2)
            φ .= max.(φ, tanh.((radius .- distance) / (2 * system.epsilon)))
        end
        
    else
        throw(ArgumentError("Unknown initial type: $initial_type"))
    end
    
    return φ, X, Y
end

"""
    create_laplacian_operator(system::AllenCahnSystem)

Create Laplacian operator for the given boundary conditions.
"""
function create_laplacian_operator(system::AllenCahnSystem)
    nx, ny = system.grid_size
    Lx, Ly = system.domain_size
    
    if system.boundary_condition == :periodic
        # Use FFT for periodic boundaries
        kx = 2π * fftfreq(nx, nx/Lx)
        ky = 2π * fftfreq(ny, ny/Ly)
        Kx = kx' .* ones(ny)
        Ky = ones(nx)' .* ky
        laplacian_k = -(Kx.^2 + Ky.^2)
        return laplacian_k, :spectral
        
    else
        # Use finite differences for other boundary conditions
        dx, dy = Lx/nx, Ly/ny
        return (dx, dy), :finite_difference
    end
end

"""
    apply_laplacian!(φ, ∇²φ, laplacian_op, method, system::AllenCahnSystem)

Apply Laplacian operator to the phase field.
"""
function apply_laplacian!(φ, ∇²φ, laplacian_op, method, system::AllenCahnSystem)
    if method == :spectral
        φ_fft = fft(φ)
        ∇²φ_fft = laplacian_op .* φ_fft
        ∇²φ .= real.(ifft(∇²φ_fft))
        
    elseif method == :finite_difference
        dx, dy = laplacian_op
        nx, ny = size(φ)
        
        if system.boundary_condition == :neumann
            # Neumann boundary conditions: zero flux
            for j in 1:ny, i in 1:nx
                i_prev = max(1, i-1)
                i_next = min(nx, i+1)
                j_prev = max(1, j-1)
                j_next = min(ny, j+1)
                
                # Handle boundaries
                if i == 1
                    i_prev = 2  # Neumann condition
                elseif i == nx
                    i_next = nx - 1
                end
                
                if j == 1
                    j_prev = 2
                elseif j == ny
                    j_next = ny - 1
                end
                
                ∇²φ[i,j] = (φ[i_prev,j] - 2*φ[i,j] + φ[i_next,j]) / dx^2 + 
                           (φ[i,j_prev] - 2*φ[i,j] + φ[i,j_next]) / dy^2
            end
        end
    end
end

"""
    evolve_allen_cahn!(φ, φ_new, ∇²φ, laplacian_op, method, system::AllenCahnSystem, dt)

Perform one time step evolution of the Allen-Cahn equation.
"""
function evolve_allen_cahn!(φ, φ_new, ∇²φ, laplacian_op, method, system::AllenCahnSystem, dt)
    # Compute Laplacian
    apply_laplacian!(φ, ∇²φ, laplacian_op, method, system)
    
    # Compute potential derivative
    f_prime = potential_derivative(φ, system)
    
    # Allen-Cahn evolution: ∂φ/∂t = M[ε²∇²φ - f'(φ)]
    φ_new .= φ .+ dt * system.mobility * (system.epsilon^2 * ∇²φ .- f_prime)
end

"""
    calculate_allen_cahn_energy(φ, ∇²φ, laplacian_op, method, system::AllenCahnSystem)

Calculate the total free energy of the Allen-Cahn system.
"""
function calculate_allen_cahn_energy(φ, ∇²φ, laplacian_op, method, system::AllenCahnSystem)
    nx, ny = system.grid_size
    dx = system.domain_size[1] / nx
    dy = system.domain_size[2] / ny
    
    # Compute gradient energy: ε²/2 |∇φ|²
    if method == :spectral
        φ_fft = fft(φ)
        kx = 2π * fftfreq(nx, nx/system.domain_size[1])
        ky = 2π * fftfreq(ny, ny/system.domain_size[2])
        Kx = kx' .* ones(ny)
        Ky = ones(nx)' .* ky
        
        gradient_energy_k = 0.5 * system.epsilon^2 * (Kx.^2 + Ky.^2) .* abs.(φ_fft).^2
        gradient_energy = real(sum(gradient_energy_k)) * dx * dy / (nx * ny)
        
    else
        # Finite difference approximation
        ∇φ_x, ∇φ_y = finite_difference_gradient(φ, dx, dy)
        gradient_energy = 0.5 * system.epsilon^2 * sum(∇φ_x.^2 + ∇φ_y.^2) * dx * dy
    end
    
    # Potential energy
    potential_energy_density = potential_energy(φ, system)
    potential_total = sum(potential_energy_density) * dx * dy
    
    total_energy = gradient_energy + potential_total
    return total_energy, gradient_energy, potential_total
end

"""
    finite_difference_gradient(φ, dx, dy)

Compute gradient using finite differences.
"""
function finite_difference_gradient(φ, dx, dy)
    nx, ny = size(φ)
    ∇φ_x = zeros(nx, ny)
    ∇φ_y = zeros(nx, ny)
    
    # Central differences with periodic boundaries
    for j in 1:ny, i in 1:nx
        i_next = i == nx ? 1 : i + 1
        i_prev = i == 1 ? nx : i - 1
        j_next = j == ny ? 1 : j + 1
        j_prev = j == 1 ? ny : j - 1
        
        ∇φ_x[i,j] = (φ[i_next,j] - φ[i_prev,j]) / (2*dx)
        ∇φ_y[i,j] = (φ[i,j_next] - φ[i,j_prev]) / (2*dy)
    end
    
    return ∇φ_x, ∇φ_y
end

"""
    run_allen_cahn(system::AllenCahnSystem; t_final=1.0, dt=nothing, 
                   initial_type=:random, save_interval=0.1, verbose=true)

Run Allen-Cahn equation simulation.
"""
function run_allen_cahn(system::AllenCahnSystem; t_final=1.0, dt=nothing,
                        initial_type=:random, save_interval=0.1, verbose=true)
    
    if verbose
        println("⚗️  Allen-Cahn Phase Field Simulation")
        println("====================================")
        println("Mobility M = $(system.mobility)")
        println("Interface width ε = $(system.epsilon)")
        println("Potential: $(system.potential)")
        println("Domain: $(system.domain_size)")
        println("Grid: $(system.grid_size)")
        println("Boundary: $(system.boundary_condition)")
        println()
    end
    
    # Initialize system
    φ, X, Y = initialize_allen_cahn(system, initial_type)
    nx, ny = system.grid_size
    
    # Create differential operators
    laplacian_op, method = create_laplacian_operator(system)
    
    # Determine stable time step
    if dt === nothing
        if method == :spectral
            max_k2 = maximum(abs.(laplacian_op))
            dt = 0.1 / (system.mobility * system.epsilon^2 * max_k2)
        else
            dx, dy = laplacian_op
            dt = min(dx^2, dy^2) / (4 * system.mobility * system.epsilon^2) * 0.8
        end
        dt = min(dt, 0.01)  # Upper bound
    end
    
    n_steps = Int(t_final / dt)
    n_save = max(1, Int(save_interval / dt))
    
    if verbose
        println("Time integration:")
        println("  dt = $(round(dt, digits=6))")
        println("  Steps: $n_steps")
        println("  Save interval: $(save_interval)")
        println()
    end
    
    # Storage
    times = Float64[]
    snapshots = Array{Float64,3}(undef, nx, ny, 0)
    energy_history = Float64[]
    interface_length_history = Float64[]
    
    # Work arrays
    φ_new = similar(φ)
    ∇²φ = similar(φ)
    
    # Save initial condition
    t = 0.0
    push!(times, t)
    snapshots = cat(snapshots, reshape(φ, nx, ny, 1), dims=3)
    
    # Calculate initial energy
    total_energy, grad_energy, pot_energy = calculate_allen_cahn_energy(φ, ∇²φ, laplacian_op, method, system)
    push!(energy_history, total_energy)
    
    # Calculate initial interface length
    interface_length = estimate_interface_length(φ, system)
    push!(interface_length_history, interface_length)
    
    if verbose
        println("Initial conditions:")
        println("  Total energy: $(round(total_energy, digits=6))")
        println("  Interface length: $(round(interface_length, digits=4))")
        print("Progress: ")
    end
    
    # Time evolution
    for step in 1:n_steps
        # Time step
        evolve_allen_cahn!(φ, φ_new, ∇²φ, laplacian_op, method, system, dt)
        φ, φ_new = φ_new, φ
        t += dt
        
        # Save data
        if step % n_save == 0
            push!(times, t)
            snapshots = cat(snapshots, reshape(φ, nx, ny, 1), dims=3)
            
            # Calculate diagnostics
            total_energy, grad_energy, pot_energy = calculate_allen_cahn_energy(φ, ∇²φ, laplacian_op, method, system)
            push!(energy_history, total_energy)
            
            interface_length = estimate_interface_length(φ, system)
            push!(interface_length_history, interface_length)
            
            if verbose
                progress = step / n_steps * 100
                print("$(Int(round(progress)))% ")
            end
        end
    end
    
    if verbose
        println("\n✅ Simulation completed!")
        final_energy = energy_history[end]
        energy_change = final_energy - energy_history[1]
        println("Final energy: $(round(final_energy, digits=6))")
        println("Energy change: $(round(energy_change, digits=6)) (should be ≤ 0)")
        println()
    end
    
    return (
        times = times,
        φ_final = φ,
        snapshots = snapshots,
        energy_history = energy_history,
        interface_length_history = interface_length_history,
        X = X,
        Y = Y,
        system = system,
        dt = dt
    )
end

"""
    estimate_interface_length(φ, system::AllenCahnSystem)

Estimate the total interface length using the gradient method.
"""
function estimate_interface_length(φ, system::AllenCahnSystem)
    nx, ny = size(φ)
    dx = system.domain_size[1] / nx
    dy = system.domain_size[2] / ny
    
    # Compute gradient magnitude
    ∇φ_x, ∇φ_y = finite_difference_gradient(φ, dx, dy)
    ∇φ_magnitude = sqrt.(∇φ_x.^2 + ∇φ_y.^2)
    
    # Interface length ≈ ∫ |∇φ| dx dy
    interface_length = sum(∇φ_magnitude) * dx * dy
    return interface_length
end

"""
    analyze_interfaces(results; verbose=true)

Analyze interface properties and dynamics in Allen-Cahn simulation results.
"""
function analyze_interfaces(results; verbose=true)
    if verbose
        println("🔬 Allen-Cahn Interface Analysis")
        println("================================")
    end
    
    φ_final = results.φ_final
    system = results.system
    times = results.times
    energy_history = results.energy_history
    
    # Basic statistics
    φ_mean = mean(φ_final)
    φ_std = std(φ_final)
    φ_min, φ_max = extrema(φ_final)
    
    if verbose
        println("Phase field statistics:")
        println("  Mean: $(round(φ_mean, digits=4))")
        println("  Std:  $(round(φ_std, digits=4))")
        println("  Range: [$(round(φ_min, digits=4)), $(round(φ_max, digits=4))]")
    end
    
    # Energy analysis
    initial_energy = energy_history[1]
    final_energy = energy_history[end]
    energy_dissipation = initial_energy - final_energy
    
    # Energy dissipation rate
    if length(times) > 1
        energy_rate = (energy_history[end] - energy_history[end-1]) / (times[end] - times[end-1])
    else
        energy_rate = 0.0
    end
    
    if verbose
        println("Energy analysis:")
        println("  Initial energy: $(round(initial_energy, digits=6))")
        println("  Final energy: $(round(final_energy, digits=6))")
        println("  Energy dissipated: $(round(energy_dissipation, digits=6))")
        println("  Final dissipation rate: $(round(energy_rate, digits=8))")
    end
    
    # Interface analysis
    final_interface_length = estimate_interface_length(φ_final, system)
    interface_evolution = results.interface_length_history
    
    if length(interface_evolution) > 1
        interface_change = interface_evolution[end] - interface_evolution[1]
        interface_rate = (interface_evolution[end] - interface_evolution[end-1]) / (times[end] - times[end-1])
    else
        interface_change = 0.0
        interface_rate = 0.0
    end
    
    if verbose
        println("Interface analysis:")
        println("  Final interface length: $(round(final_interface_length, digits=4))")
        println("  Interface change: $(round(interface_change, digits=4))")
        println("  Final interface rate: $(round(interface_rate, digits=6))")
    end
    
    # Phase separation analysis
    # Count connected components (simplified)
    positive_fraction = count(φ_final .> 0) / length(φ_final)
    negative_fraction = count(φ_final .< 0) / length(φ_final)
    interface_fraction = count(abs.(φ_final) .< 0.1) / length(φ_final)
    
    if verbose
        println("Phase analysis:")
        println("  Positive phase: $(round(positive_fraction*100, digits=1))%")
        println("  Negative phase: $(round(negative_fraction*100, digits=1))%")
        println("  Interface region: $(round(interface_fraction*100, digits=1))%")
        println()
    end
    
    return (
        statistics = (mean=φ_mean, std=φ_std, min=φ_min, max=φ_max),
        energy = (initial=initial_energy, final=final_energy, 
                 dissipated=energy_dissipation, rate=energy_rate),
        interface = (length=final_interface_length, change=interface_change, 
                    rate=interface_rate),
        phases = (positive=positive_fraction, negative=negative_fraction,
                 interface=interface_fraction)
    )
end

"""
    track_energy_dissipation(results; plot_results=true)

Create detailed energy dissipation analysis.
"""
function track_energy_dissipation(results; plot_results=true)
    times = results.times
    energy_history = results.energy_history
    interface_history = results.interface_length_history
    
    println("📈 Energy Dissipation Analysis")
    println("=============================")
    
    # Calculate energy dissipation rate
    if length(times) > 2
        dt_avg = (times[end] - times[1]) / (length(times) - 1)
        energy_rates = diff(energy_history) / dt_avg
        
        println("Energy dissipation statistics:")
        println("  Total dissipated: $(round(energy_history[1] - energy_history[end], digits=6))")
        println("  Average rate: $(round(mean(energy_rates), digits=8))")
        println("  Final rate: $(round(energy_rates[end], digits=8))")
    end
    
    # Interface evolution
    if length(interface_history) > 2
        interface_rates = diff(interface_history) / dt_avg
        println("Interface evolution:")
        println("  Total change: $(round(interface_history[end] - interface_history[1], digits=4))")
        println("  Average rate: $(round(mean(interface_rates), digits=6))")
    end
    
    if plot_results && length(times) > 2
        # Create energy plot
        p1 = plot(times, energy_history, xlabel="Time", ylabel="Energy", 
                 title="Energy Evolution", linewidth=2, legend=false)
        
        p2 = plot(times, interface_history, xlabel="Time", ylabel="Interface Length",
                 title="Interface Evolution", linewidth=2, legend=false, color=:red)
        
        # Combine plots
        energy_plot = plot(p1, p2, layout=(2,1), size=(600, 500))
        
        return energy_plot
    end
    
    return nothing
end

# Demonstration and validation examples
if abspath(PROGRAM_FILE) == @__FILE__
    println("⚗️  Allen-Cahn Equation - Research Examples")
    println("=" * 50)
    
    # Example 1: Interface motion
    println("\n1. Interface Motion (Mean Curvature Flow)")
    system1 = AllenCahnSystem(mobility=1.0, epsilon=0.05)
    results1 = run_allen_cahn(system1, t_final=0.5, 
                             initial_type=:circular_interface)
    analysis1 = analyze_interfaces(results1)
    energy_plot1 = track_energy_dissipation(results1)
    
    # Example 2: Phase separation
    println("\n2. Phase Separation Dynamics")
    system2 = AllenCahnSystem(mobility=1.0, epsilon=0.02,
                             domain_size=(4.0, 4.0), grid_size=(256, 256))
    results2 = run_allen_cahn(system2, t_final=2.0,
                             initial_type=:random)
    analysis2 = analyze_interfaces(results2)
    
    # Example 3: Multiple droplets coarsening
    println("\n3. Droplet Coarsening")
    system3 = AllenCahnSystem(mobility=0.5, epsilon=0.03,
                             domain_size=(3.0, 3.0), grid_size=(192, 192))
    results3 = run_allen_cahn(system3, t_final=3.0,
                             initial_type=:droplets)
    analysis3 = analyze_interfaces(results3)
    
    println("\n🎉 All Allen-Cahn examples completed!")
    println("Validated energy dissipation and interface dynamics.")
end